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Abstract 

Two-state models provide phenomenological descriptions of many different systems, ranging from 
physics to chemistry and biology. We investigate work fluctuations in an ensemble of two-state systems 
driven out of equilibrium under the action of an external perturbation. We calculate the probability den- 
sity P^iyV) that a work equal to W is exerted upon the system (of size N) along a given non-equilibrium 
trajectory and introduce a trajectory thermodynamics formalism to quantify work fluctuations in the 
large- limit. We then define a trajectory entropy SMiW) that counts the number of non-equilibrium 
trajectories P^iW) = eyi^^SNiW) /^bT) with work equal to W and characterizes fluctuations of work 
trajectories around the most probable value T^™p. A trajectory free- energy .F/v(VF) can also be defined, 
which has a minimum sX W = W\ this being the value of the work that has to be efficiently sampled 
to quantitatively test the Jarzynski equality. Within this formalism a Lagrange multiplier is also intro- 
duced, the inverse of which plays the role of a trajectory temperature. Our general solution for PjsiiW) 
exactly satisfies the fluctuation theorem by Crooks and allows us to investigate heat-fluctuations for a 
protocol that is invariant under time reversal. The heat distribution is then characterized by a Gaussian 
component (describing small and frequent heat exchange events) and exponential tails (describing the 
statistics of large deviations and rare events). For the latter, the width of the exponential tails is related 
to the aforementioned trajectory temperature. Finite-size effects to the large-A theory and the recovery 
of work distributions for flnite N are also discussed. Finally, we pay particular attention to the case of 
magnetic nanoparticle systems under the action of a magnetic field H where work and heat fiuctuations 
are predicted to be observable in ramping experiments in micro-SQUIDs. 
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1 Introduction 



There has been recent interest in the experimental measure of work fluctuations and the test of the 
so-caUed fluctuation theorems. Initially proposed in the context of sheared systems in a steady state P 
several versions of such theorems have been derived [2]. In particular, speciflc identities have been 
obtained in the context of stochastic systems that show how it is possible to recover the equilibrium 
free-energy change in a reversible transformation by exponential averaging over many non-equilibrium 
trajectories that start at equilibrium jSHHISl- Let us consider a system initially in equilibrium in contact 
with a thermal bath (at temperature T) that is submitted to an isothermal perturbation according to 
a given protocol. Work fluctuations (WF) refer to the fact that the work W exerted upon the system 
depends on the particular non-equilibrium trajectory followed by the system. As the initial conflgura- 
tion or the trajectory are stochastic, the value of the work W changes among different trajectories, all 
generated with the same perturbation protocol. Transient violations (TV) of the second law refer to the 
fact that, among all possible WF, a fraction of them absorb heat from the bath that is transformed into 
work. Taken individually, these rare trajectories violate the Clausius inequality, Q < TAS, where Q is 
the heat supplied from the bath to the system and AS* is the change in the entropy, a state function 
defined through the transformation. In a transformation cycle AS" = these TV satisfy Q > i.e., they 
can absorb a net amount of heat from the bath during the cycle. In terms of the dissipated work W^is, 
the Clausius relation can be expressed in the following form, 

W^dis = W- Wre. = TAS -Q>0 (1) 

In this expression W is the total work exerted upon the system. According to the first law of thermody- 
namics (conservation of the energy) W is given hy W = AE — Q where AE is the change in the internal 
energy, W^ev is the reversible work (identical to the free-energy change AF = AE — TAS). Both the heat 
Q and H^dis (or W) are trajectory dependent, however AS" and Wrev are both trajectory independent 
as they are state functions, only dependent on the initial and final states. The Clausius inequality 
has to be understood as a result that is valid after averaging the fiuctuating quantities Q and W over 
an infinite number of trajectories (in what follows we will denote this average by (..)). The second law 
reads Wdis ^ and TV of the second law refer to the existence of trajectories where PFdis < 0. From this 
point of view, TV are just WF characterized by the fact that W^is < 0. The interest in studying TV is 
that these describe large deviations of the work that have to be sampled in order to recover equilibrium 
free-energy differences from non-equilibrium measurements [HI. 

The steadily increasing development of nanotechnologies during the last decade has made WF experi- 
mentally accessible. Recent experiments on single RNA hairpins unfolded under the action of mechanical 
force [7] and micro-sized beads trapped by laser tweezers and moved through a solvent |S1 have provided 
a first quantitative estimate of WF and TV. Related measurements include the experimental verifica- 
tion of the Gallavotti-Cohen fluctuation theorem in Rayleigh-Bernard convection [0] and turbulent flows 
|lUj . This research is potentially very interesting as it leads to new insights about the physical processes 
occurring at the nanoscale, a frontier that marks the onset of complex organization of matter A 
characteristic of WF is that they are quickly suppressed as the system size or the time window of the 
measurement increase. 

The central quantity describing WF is the work probability distribution Pn{W) {N stands for the 
system size), PN{W)dW being the fraction of non-equilibrium trajectories with work between W and 
W + dW. The knowledge of this quantity is important for what it tells us about the mathematical 
form of the tails of the distribution, relevant to understand the importance of large deviations of work 
values respect to the average value. A precise knowledge of the form of the tails in that distribution 
gives us hints about how many experiments need to be done in order to recover equilibrium quantities 
from non-equilibrium experiments. In this work we investigate an ensemble of two-level systems as an 
explicit example where P]\f{W) can be analytically computed in the large- approach using a path 
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integral method. This approach allows us to exactly derive several exact results describing work and 
heat fluctuations in the system in the large-A^ limit but also for finite N. The most important result 
in the paper is the introduction of a trajectory thermodynamics formalism, the key quantity being the 
trajectory entropy Sn{W). This allows us to infer several quantities such as the trajectory free-energy 
^n(W) and the trajectory temperature \{W), the latter being a Lagrange multiplier that plays the role 
of the inverse of a temperature, an intensive variable related to the statistics of large deviations or tails in 
the work and heat distributions. Two-state models represent a broad category of systems where WF and 
TV can be predicted to be experimentally observable making the present calculations relevant as they 
might allow a detailed comparison between theory and experiments. In particular, we propose magnetic 
nanoparticles as excellent candidate systems to experimentally test the present theory. 

The plan of the paper is as follows. In Sec. [21 we describe the model and the large- approach. 
In Sec. 01 we develop the trajectory thermodynamics formalism that allows us to reconstruct the work 
distribution, define a trajectory entropy S]\r{W) = log(PAr(VF)) and a trajectory free-energy J^j\[(W). 
In Sec. m we show how the saddle-point equations derived in Sec. |21 can be numerically solved. The 
dependence of the main parameters of the theory (most probable work VF™^, transient violations work 
and fluctuation-dissipation ratio R) on the field protocol are discussed in Sec. 14.11 Within the 
formalism it is then possible to show, Sec. |SJ that the entropy per particle s{w) {w being the work W 
per particle) exactly satisfies the fluctuation theorem by Crooks. Moreover, it is possible to infer the 
shape of the tails in the work distribution from the sole knowledge of the Lagrange multiplier conjugated 
to the trajectory entropy, X{w), that plays the role of the inverse of a temperature (what we call the 
trajectory temperature) in the formalism. In Sec. IHl we study heat fluctuations in the model. We show 
the existence of two sectors in the heat distribution that are described by a Gaussian central part 
(corresponding to small and most probable deviations) and two exponential tails (corresponding to large 
and rare deviations) showing the presence of intermittent heat fluctuations in the theory. In Sec. [7| we 
discuss finite-size corrections to the large- theory and how Pn{W) for finite can be reconstructed 
using the results from the large- approach. Particular emphasis is finally placed in Sec.|Hlin the case of 
magnetic nanoparticle systems where WF are predicted to be experimentally observable and described 
by the present theory. Sec. IHl presents the conclusions. 



2 Ensemble of two-state systems: the large- approach 

A broad category of systems can be modeled by an ensemble or collection of independent two-state sys- 
tems. These offer realistic descriptions of electronic and optical devices that can function in two different 
configurations, atoms in their ground and excited states, magnetic particles whose magnetic moment can 
point in two directions, or biomolecules in their native and unfolded states, among others. Throughout 
the paper, and in view of the possible experimental implications, we will adopt the nomenclature of 
magnetic systems. A particle i in the ensemble (1 < « < A^) has magnetic moment /i and can point 
in two directions according to the sign of the spin cij = ±1. A given configuration in the ensemble is 
specified by a string of spin values C = {ai, 02, ■■, cr^}. In the presence of an external field H, the energy 
of a configuration C is given by 

TV 

E{C) = -iiHM{C) = -i-iHY^Ui , (2) 

i=l 

= X^i^i "^i being the total magnetization of the system. The transition rates for individual particles 
will be denoted as p^^{H),p'^"^'^{H) to indicate the transitions 0" = — 1— >(t' = 1 and a = 1 ^ a' = —1 
respectively. These rates satisfy detailed balance, therefore p^^{H) /p^°'"^{H) = exp{—2Pij,H) where 
P = l/ksT, T being the bath temperature and ks the Boltzmann constant. The overall transition rate 
is given by p^°^{H) = p^^{H) + p'^"^'^{H). Although it is possible to introduce structural disorder in the 
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ensemble (e.g. by allowing /u or p^^{H) to be a random quenched variable), in the following analysis we 
will restrict us to the non-disordered or mono-disperse case. 

Let the system be prepared at t = in an equilibrium state at an initial value of the field Hq = Hi 
and let us consider an external isothermal perturbation that changes the field H according to a protocol 
function H{t). Throughout this paper we will denote this non-equilibrium process a ramping experiment. 
If the variation is slow enough then the process is quasi-static and the system goes through a sequence 
of equilibrium states. However, if the rate H is large compared to the relaxation time of the particle 
then the magnetization M = J2iLiC!'i does not follow the equilibrium curve Meq{H) = N ta,nh(P fiH) . 
To specify a trajectory it is then convenient to discretize time in A^^ time-steps of duration At each and 
take the continuous-time limit At ^ 0, A^^ ^ oo (with the total time t = NgAt fixed) at the end. The 
perturbation protocol is specified by the sequence of values {Hk;l < k < Ns} and a trajectory T is 
defined by the sequence of configurations T = {Ck; I < k < Ng) where Ck = {erf; I < i < N} is the 
configuration at time t = kAt. The total work exerted upon the system along a given trajectory is given 
by 0, 

Ns-l 

W{T) = -l^J2 Mk+iiHk+i - Hk) (3) 

Mk = J2iLi being the magnetization at time-step k. The dissipated work for a given trajectory is 
the difference between the total work and the reversible one, W^is = W — Wrev where VFrev = AF is 
the change in equilibrium free-energy between the initial and final values of the field. The free energy 
is given by F{H) = —NksT log{2cos]i{f3 fj,H)). To quantify WF we have to compute the probability 
distribution for the total work measured over all possible non-equilibrium trajectories, 

Ns-l 

PNiW) = J2p{T)6{W - W{T)) = J2 p{T)SiW + Mk+iiHk+i - Hk)) , (4) 

T {^fc} fc=0 

where p{T) denotes the probability of a given trajectory. The subindex N in Pfyf{W) is written to em- 
phasize the dependence of the distribution on the size of the system. Pn{W) is computed using the Bayes 
formula p{T) = Ylk=o^ 1kiWi~^^}\Wi})PoiWi})^ where QkiWjlW}) denotes the transition probability 
to go from {fj} to {a'} at time-step k, and po({o'i'}) is the initially equilibrated (i.e. Boltzmann-Gibbs) 
distribution. Evaluation of the integral Q) requires the following steps: 1) trace out spins in the sum; 
2) insert the factorized expression for p(T); 3) use the integral representation for the delta function 
6{x) = (l/27r) dx exp(zAx) and 4) insert the following factor, 

Ns-l I noo N 

l=n^/ d-fkdMkexp{ijkiMk-J2a'^)) . (5) 
k=o ^'^J-oo 

After some manipulations this leads to the following expression for the work probability distribution (up 
to some unimportant multiplicative terms), 

„ Ns-l 

PAr(l^) oc d\ Yl {d'ykdmk)ew{A{w,\{7k},{'mk})) (6) 
fc=o 

where A is the saddle-point function, w = W/N,mk = Mk/N (throughout the paper we will use small 
case letters to refer to intensive quantities). The function a = A/N is given by, 

Ns-l Ns 

a{w, A, {7fc}, {ruk}) = -X(w + ^ mk+i{Hk+i - Hk)^ - Ikmk + 

fc=o fc=o 

^ \og{uk+i) + \og{vk+i)) + logie^'p-^m + e-^Y'-^^Hi)) . (7) 

fc=0 
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The terms Uk,Vk are given by, 



Uk+i = exp(7fe+i)(l - pt°n + exp(-7fc+i)pf-'^ (8) 
Vk+i = exp(7fc+i)p^P + exp(-7fc+i)(l - p^)- (9) 

with the boundary condition 'j^^ = 0. The quantities p^^{Hi)^p^'^{Hi) are the transition rates at time 
s = 0, and we are assuming that at the initial condition the system is in thermal equilibrium. In the 
continuous-time limit © becomes a path integral over the variable A and the functions ^{t),m{t) with, 

t ^ 

+ 



where 



a{w, X,j{s), 171(3)) = —X(w + fi J 'm{s)H{s)ds^ 
i [m{s){2j{s) + c{s)) + d{s))ds + log(e^Wp"P(/fi) + e"^(°)/°"'X^i)) (10) 



c(,) =/°--(s)(exp(-27(.)) - 1) -p"P(.)(exp(27(.)) - 1) (11) 
d{s) = /°-'^(s)(exp(-27(s)) - 1) +p-P(.)(exp(27(.)) - 1) . (12) 

As we are interested in the crossover to the large- regime we can estimate the integral © by using the 
saddle-point method. For each value of the work trajectory w the dominant contribution is given by the 
solution of the functional equations, 



w + fi f m{s)H{s)ds = (13) 
Jo 



5a '■^ 



m{s) + m{s)p'°'{s) - (p"P(s) - /°""^(s)) + m{s)d{s) + c{s) = (14) 



5^(3 



7(^) - ^^^His) + Us) = (15) 



6m{s) ' ' ' 2 

with the boundary conditions 

7(t) = ; m(0) = tanh(7(0) + pfiHi) . (16) 

Note that the boundary conditions are a bit special as causality is broken. The function 7(5) has the 
boundary condition located at the final time s = t while the boundary condition for m{s) is located 
at the initial time s = 0. These equations can be numerically solved in general and analytically solved 
only partially and for some particular cases (e.g. in the case where the rate H is constant). Before 
presenting detailed numerical solutions to these equations we should point out several general aspects 
of such solutions. At first we note how, for a given value of A, Eq. (|15|) together with the boundary 
condition 7(t) = can be solved giving the solution 7a(-s), the subindex A emphasizing the dependence 
of this solution on the parameter A. Inserting this result in (|14|) and using the boundary condition (|l()j) 
we get the solution mx{s). Finally, insertion of mx{s) in ()13() gives a value for the work w{\). This last 
relation can then be inverted to give X^w) and from it, the solutions "yx{s),mx{s) will also depend on 
the value of w. To better emphasize this dependence we will denote by X{w),^w{s),mw{s) the solutions 
of (|13I14I15|) for a given value of w and 

s{w) = a{w,X{w),-/w{s),m^^{s)) (17) 

the corresponding extremal value of a. We will also make explicit the lo-dependence in the time-dependent 
quantities c(s), d{s) in (|11I12() and denote them by Cw{s), dw{s) respectively. Furthermore, we can define 
the trajectory entropy Sn(W), 

Pn{W) = eMSNiW)) . (18) 
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In the large- A/" limit, from (|f)|17)l we have 

s{w) = lim ^^[^'> with W = Nw , (19) 

the function s{w) playing the role of a trajectory entropy per particle that counts the density of tra- 
jectories per particle with work equal to w. This means that, for A'" finite, ^N{w)dw = exp{Ns{w))d'w 
is approximately proportional to the fraction of trajectories with work between w and w + dw. From 
(|18I19)) an approximate expression for the work probability distribution can be written, 

p ( \^ '^n{w) ^ exp{Ns{w)) 
^ LT'''^N{w')dw' r^^''exp{Ns{w'))dw' 



where Wmm and Wmax are the minimum and maximum possible values of the work. Clearly, from (jSJ these 
values are given by tUmax = — 'U^min = t^-iH f — Hi) where Hj = H{t) is the final value of the magnetic field. 
The subindex in Pn{w) and ^n{w) emphasize the dependence of these quantities on the size of the 
system. Finally, we note that, albeit the solutions (|13I14I15|) have been obtained using the saddle-point 
approximation (only valid for large N) the final result (|2()|) can be very accurate for small values of A'^. 
This result, that at first glance may appear striking, is just consequence of the non- interacting character 
of the Hamiltonian (jj}. This point is discussed in more detail in Sec. [3 There we show that, albeit (|2nj) 
is only approximate for finite A^, the cumulants that we can extract from s{w) are exact for any A^. This 
allows us to exactly reconstruct the finite A^ distribution from the sole knowledge of s{w). 

The action A = Na in ()1U() could be used (employing Monte Carlo algorithms) to generate trajectories 
according to their probability Pn{w) ^- Inserting (|T5|) in (fTTHl we get, 

s{w) = —\w + - / du,{s)ds . (21) 
2 Jo 

The value w for which s{w) is maximum yields the most probable work {w = w^"^) among all trajectories. 
This can be evaluated using the equation 

s'{w) = — = -X{w) , (22) 

where we have used the chain rule together with the extremum conditions pHll4ll5j) as well as We 
will see later in SeclHlthat the Lagrange multiplier X{w) is related to the inverse of a new energy scale or 
temperature that describes the tails of the work distribution. This quantity is of much current interest 
as it describes the statistics of rare events and large deviations of work values from the average which 
are observable in small systems. The extremum solution of (|22|) can then be written as A(u;™p) = 0, 



ds{w) 



dw 



= or A(u;°>p) = . (23) 

UI=UI™P 



This solution solves (|13I14I15|) giving 7u;„jp(s) = c^mp(s) = dw^p{s) = 0. Eas. H13ll4|) then give the 
solution for the most probable trajectory (usually derived using standard statistical methods), 

mis) = -m(s)p*°*(s) + (p"P(s) - p'^°^'^(s)) . (24) 

The reversible process is a special case (only valid for slow enough perturbation protocols) and corre- 
sponds to m(s) = or 

m{s) = (p"P(s) - /°"''(s))/p*°*(s) = tanhip fiH{s)) . (25) 



"'^The easiest procedure then would be to start from an initial trajectory j(s),m{s) (satisfying the boundary conditions 
to(0) — tanh(7(0) + PnHi);-f{t) — 0) and perform successive "local" updates along the trajectory and accepting the moves 
according to the change in the action A (by using an algorithm that satisfies detailed balance, as defined by the action A, and 
respects the boundary conditions). 
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3 Trajectory thermodynamics formalism 



Prom the trajectory entropy s{w) we can construct a trajectory free-energy ^(w) useful to predict under 
which conditions TV are properly sampled and fluctuation theorems can be quantitatively verified. For 
this we consider the Jarzynski equality j^ij 



that we can write as, 

where we used (|18j) and we have defined the trajectory free-energy, 

J^NiW) = W - keTSNiW) . (28) 
In the large- limit, using (|19j) . we can write 



dw exp|^- 



where 

j:{w) = w- kBTs{w) (30) 

is a trajectory free-energy (per particle) that depends on the particular value of the work w. Evaluating 
the integral by the steepest descent method and using H21() we obtain the thermodynamic relations. 



ds{w) 



ksT dw 

ksT 



= -\{w^) (31) 



kuT r 

J^{w^) = AF/N = wre^ = w^ - kBTs{w^) = — ^ / d^t{s)ds (32) 

2 JO 

Using the definition (|30|) together with (|23I3H) we have the relations. 



dw 
d!F{w 



= 1 (33) 
, = (34) 



dw 

i.e. the entropy has a maximum at it; = w™^ and the free energy has a minimum ok, w = w'^ . These 
relations bear similarity to those considered in thermodynamics but now applied to work trajectory 
values. For the case of the canonical ensemble the quantities s{w),J^{w),w play the role of the standard 
entropy, free energy and internal energy while X{w) is the intensive variable corresponding to the inverse 
of a temperature. 

A graphical construction of the relations ()31I32() is shown in Figure ^ This figure illustrates how the 
most important quantities Wrev^w'^^'^jW^w are related to each other. In particular, w = J dwwPj\[{w) is 
expected to differ from albeit that difference can be small for highly symmetric distributions. 

The difference between w^ and ?i)™P indicates that the average (|29|1 is properly weighed whenever 
trajectories with work values around w'^ are sampled. This result indicates that proper sampling of non- 
equilibrium work values around w^ is required to derive equilibrium free-energies from non-equilibrium 
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Figure 1: Diagrams showing the different relevant quantities in the trajectory thermodynamics formalism. Upper 
panel: Trajectory entropy s{w) related by (|2flj) to the density of trajectories with work equal to w. Middle panel: 
Trajectory free-energy J-{w) = w — kBTs{w). Lower panel: Lagrange multiplier X{w). Six are the most relevant work 
values: Wmax and ifmin for the maximum and minimum values of the work; w"^'^ the most probable work value given 
by s'(t(;™P) = A(?i;™P) = or !F'{w'^^) = 1; the value of the work that has to be sampled to recover free energies 
from non-equilibrium work values using the Jarzynski equality (|26|) . This is given by s'{w^) = —X{w^) = l/ksT or 
^'{w'^) = 0; Wrev = F{Hf) — F{Hi) the reversible work; and VMi^ = J^^^^i'^ ~ Wiev)PN{w)dw the average dissipated 
work. They are related by Wmm < < w'^^ < tfmax while the second law of thermodynamics imposes w^is > 0. 
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measurements by using the Jarzynski equality. A proper sampling of work values around can be 
guaranteed when, out of the total number of trajectories, a finite fraction of work trajectory values in 
the vicinity of is observed. From a practical point of view this means that the histogram of work 
values must extend down to . If this is not achieved, then the exponential average performed over a 
finite number of non-equilibrium experiments has a bias that can be estimated in some cases jl41 115| . 
Eas. H31l32() are readily solved at the Gaussian level (i.e. assuming that Pn{w) is exactly a Gaussian 
or s{w) a quadratic function) giving = tt;™^ — a'^/ksT (u^ being the variance of the Gaussian 
work distribution). For quasi-reversible processes in the linear response regime the fluctuation- 
dissipation theorem implies ~ 2kBTw'^^ giving ~ """^di^) trajectories with negative values 
of the dissipated work that are of the order (in absolute value) of the average dissipated work must 
be sampled to quantitatively verify the validity of the Jarzynski equality. An example of such a quasi- 
reversible process, where PN{wdis) is exactly a Gaussian, is the case of a Brownian particle subjected to 
an harmonic potential and dragged in a fluid [51 IT^ [T7| . 

4 Numerical solution of the equations 

Equations (|13I14I15|) can be numerically solved in general. We assume Glauber transition rates given by 

p^^{H)=p'^\H)q{H) ■ p''°--iH)=p'°\H)il-qiH)) (35) 

with q{H) = (1 + tanh(/3/ii/))/2 and p^°'^{H) = l/T^clax{H) = a{H) corresponding to the inverse of the 
relaxation time. In this case, 

p [s) «l^i^ii2cosh(/3/xi/(s)) ^^^^ 

p [s) «i^i«ii2cosh(/3/x//(s)) ■ ^^-^ 
Inserting these expressions in (|11I12() we obtain, 

c{s) = -a(i/(.))^^^^i^MfI±^^ + aiH{s)) tanh(/3^//(.)) (38) 
cosh(p/iii (sj) 

\ /rr^ ^^COSh(27(s) +/3/ii?(s)) frrfW /QOA 

dis) = a{His)) a{His)) . (39) 

The solution of the equations consists of the following steps: 

1. Solution of 7a(s). With the boundary condition at the final time s = t, 'yx{t) = 0, Eq. (|15|) has 
to be numerically integrated backwards in time. Inserting (|36l37j) in (|15j) we obtain, 

7(s) = XfiHis) + a{H{s)) sinh(7(s))(cosh(7(s)) + sinh(7(s)) tanh(/3/ii?(s))) . (40) 

However, a direct numerical integration of this equation leads to divergences and numerical insta- 
bilities. It is then convenient to express H4U|) in terms of a new variable e(s) = l/cosh(7(s)) which 
displays smooth behavior. Equation ()4U() becomes, 

cosh(sj V 

a{H{s)) sinh(7(s))(-J-- + sinh(7(s)) tanh(/3/ii/(s)))) (41) 
e[s) 

with the boundary condition e(t) = 1. This equation can then be easily numerically integrated to 
give 7a (s) for a given value of A. 
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2. Solution of mA(s). Once the solution of (|41|) for a given value of A, "y\{s), is found, then it is 
possible to integrate ()14() to find mx{s). Because ()14() is linear its solution can be explicitly written, 

P S P S P s 

mx{s) = mx{0) exp{ Ai{u)du) + I duA2{u)eyip{ Ai{v)dv) (42) 
JO Jo Ju 

with the definitions, 

A f \ (TJ( ^^ sinh(27A(g) . . 

^^(^^ = cosh(/3^i/(.)) 

4 r \ .f./ ^^ cosh(27A(g) + 

= -a(F(.)) ^^^^^^^^^ (44) 

with the initial condition m\{Q) = tanh(7A(0) + (3^Hi). 

3. Evaluation of w, s(w), jr(w). Once ^\{s),mx{s) are known then we can evaluate w using (|13|) . 
the entropy s{w) using (|^T|) and the free energy !F{w) = w — kBTs{w). 

4. Dependence of the numerical algorithm on the sign of A. We must emphasize that the 
solution of the equations previously described only works in a sector of values of A of a given 
sign, A < 0, and leads to numerical instabilities in the other sector, A > 0, indicative that the 
transformation e(s) = 1/ cosh(7(s)) is inappropriate for A > 0. We have found a simple way out to 
this problem. It can be easily proven that the solution of (|15|) for a given value of A > is equivalent 
to the solution of that equation with the value of A with its sign reversed (—A < 0) and for the 
reversed field protocol Hr{s) = —H{s) (the subindex r stands for reversed). Eq. ((T3|) can then be 
solved and the resulting reversed solutions mr{s),^r{s),Cr{s),dr{s) give the final solutions for the 
original value of A > 0: m{s) = —mr{s),^{s) = — 7,.(s),c(s) = —Cr{s), d{s) = dr{s) (all change sign 
except d{s)). At first glance, this symmetry property might seem to be related to the content of 
the fluctuation theorem. However this relation is only apparent because the reversed process in this 
case does not correspond to the time-reversal protocol which should be instead Hr{s) = H(t — s) 
(see the discussion below in Sec. EJ- 

For the present numerical analysis, and for the sake of simplicity, we will consider a particular example 
where the ramping field H(s) changes from H{0) = Hi to H{t) = Hf at a constant rate r = H, 

H{s) = Hi + rs , r = H = ■ (45) 

We will also consider p^°^{H) = a independent of the field. This tantamount to assume that p^°^{H) 
corresponds to a microscopic attempt frequency or, rather, that the activation barrier is field independent. 
We numerically solved the equations in natural units = kgT = 1 and we have chosen a = 1 as the 
characteristic relaxation timescale of the system. Results for different values of Hi, Hf have been obtained 
by doing ramping experiments at different values of the ramping speed r. In Figure |2l we show, in the 
particular example Hi = 0,Hf = 1, the results for the magnetization trajectory solutions mx{s) and 
the Lagrange multiplier 7a (s)- These are plotted as a function of the time-dependent field H{s) for 
different values of A and for a given value of the ramping speed. In Figs. 01 lU we show several trajectory 
thermodynamics quantities at different ramping speeds (r = 0.01, 0.1, 1, 10). In the left panel of FigureOl 
we plot the magnetization for the most probable trajectories 771^=0 (s) as a function of H{s). In the 
right panel of Figure |31 and in Figure 0] we show the different trajectory thermodynamics quantities as 
a function of Wdis- the inverse temperature A(u)dis)5 the trajectory entropy s(tt'dis) and the trajectory 
free-energy J^(wdis)- 

4.1 Average and variance of the work distribution 

As has been schematically depicted in FigureEthere are different work quantities that can be of relevance 
to characterize work fluctuations. We have already defined the most probable work i(;™P and the work 
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Figure 2: Protocol with Hi = 0,Hf = 1 and ramping speed r = 1. Curves correspond to different values of 
A (A = — 5, —2, —1, —0.5, —0.2, 0., 0.2, 0.5, 1, 2, 5 from top to bottom in the upper and lower panel). Upper panel: 
Magnetization mx{s) obtained from H42|). The dashed line is the equilibrium magnetization meq{H) = tanh(ff) 
corresponding to the reversible ramping experiment r = 0. Lower panel: Lagrange multiplier 7a (s) obtained from 
(|1T|| . Note the boundary condition 7a (i) = and the presence of the most probable trajectory 7a=o(s) = 0, Vs. 




Figure 3: Protocol with Hi = 0, Hf = 1 and different ramping speeds r = 0.01, 0.1, 1, 10 (indicated by numbers along 
the continuous curves in both panels). The reversible work is w^ev = —0.433781 and Wmax = Ij'W^dis^ ~ Wmax — Wrev = 
1.433781 , uJ^j™ — Wmin ~ 'U^rev — —0.566219. Left panel: magnetization evolution for the most probable trajectories. 
The dashed line corresponds to the reversible trajectory for r ^ 0. Right panel: Lagrange multiplier A(t(;dis) for 
different ramping speeds. The intersection of the different curves with the dashed line A = gives tt;™P while the 
intersection with A = — 1 gives . The intersection of all lines at different speeds around A = —0.5 is only accidental 
(looking at a larger resolution or considering other parameters for the protocol such common crossing does not exist). 
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Figure 4: Same parameters and ramping speeds as in Figure |31 Narrower curves correspond to lower ramping 
speeds. Upper panel: Dynamical entropies s{w) plotted as functions of lUdis = w — Wrev According to the trajectory 
thermodynamics relations (|Hli:-{2|) the straight line y{wdis) = Wdis/^BT (we take fc^T = 1) is tangent to the curve 
s{wdis) at Wdis = ^^dis = — Wrev and crosses the Wdis-axis at t^dis = 0, s(?i;dis) = 0. All entropies vanish at 
"U^dis = ^dis ~ '^'^^ ~ ^rev Lower panel: Trajectory free-energy J^{W(iis)- It is identical to the equilibrium free-energy 
change AF = w^ev at w^is = ^dis- According to the same relations H31I32|) the straight line y (if dis) = "U^rev + Wdis is 
tangent to the curve J^{wdis) at Wdis = w^]^ = w"^"^ — Wrcv and crosses the tUdis-axis at Wdis = OjT^Wdis) = Wrev 
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r 



r 



Figure 5: Ramping experiments with Hi = and different values of Hj as function of the ramping speed r. Left 
panel: w^^s (continuous lines) and — tt'^jg (dashed lines) for different values of Hf indicated in the figure. The dash- 
dotted straight line corresponds to the linear response behavior Wdis = r given by (|^ . Note that zwjjjg is negative so 
we changed its sign in order to compare it with it'dis- Right panel: Fluctuation-dissipation ratio i? as a function of r 
evaluated at different values of Hf (from bottom to top, Hj = 1, 2, 3, 4, 5). 



. Another important quantity is the average work w, 

w= wPN{w)dw (46) 

where Pm{w) was defined in Q or in approximate form in (|2fl|) . In most cases (for instance, when the 
work distribution has asymmetric tails) the average work w is different from the most probable work 
w™^ . w can be lower or higher than the most probable work im^p. However, in our large-A^ theory, 
^mp _ ^j^j indistinctly both quantities in this section. We defer the discussion about 

finite-size effects in these quantities until Sec. El Another important quantity that characterizes the work 
distribution is its variance, 

=^-(lZj)2=^-(^)2 . (47) 

The average work w (or w^^) is the most relevant physical quantity that connects with classical ther- 
modynamics. The second law of thermodynamics establishes that it cannot be lower than the reversible 
work, w > w^ev However, it is clear that there can be WF such that w < w^ev These have been 
called transient violations (TV) of the second law. The relevant work value characterizing this sector 
of trajectories is given by . Clearly, W is always higher than In Figure [3 (left panel) we show 
the dependence of Wdis,w^ with the ramping speed when the field is ramped from Hi = to Hf for 
different values of Hf. It is possible to write down explicit analytic expressions for the cumulants of the 
distribution P]\f{w) in the large-A^ limit. Interestingly, and due to the non-interacting character of the 
model ((21), the cumulants derived from the large- A'^ approach are exact at all values of N, see Sec. 13 In 
particular, the first moment is given by. 



dq{H) 



2, 1 dsHis) duHiu)^l^^^^^e.^{- £ dvaiHiv))) . (48) 
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The expression for the second cumulant or variance can be obtained by expanding the function s{w) (|21|) 
up to second order with X{w) as the smaU parameter. Using the result s'{w"^^) = we get 

s{w) = + I (^^^) - w""? + 0[{w- w'^'^f] . (49) 

Prom H2U() and (|22() we obtain the relation, 



2 



1 






-1 1 dw{\) 




iV 


- dw"^ 




~ N dX 


A=0 



(50) 



We do not reproduce the details of this lengthy calculation here, the same results have been already 
obtained in a slightly different context in a previous work and in the limit of large free-energy changes 
AF as compared to ksT fT^ . 

Another interesting aspect of the present theory is that it is possible to expand the cumulants around 
the limits r ^ or r ^ oo. The former is particularly interesting because it corresponds to the so- 
called linear-response regime. In Reference |12j this regime was considered by expanding the average 
dissipated work up to linear order in the perturbation speed. By using dimensional arguments and direct 
comparison with the equivalent expression derived in the context of mechanical force jl2j we can derive 
the following result, 

where AM = M(.q{Hf) — M(,q(Hi) is the difference of equilibrium magnetizations between the initial and 
final values of the field whereas Tj-^iaxiHc = 0) = l/ptot{Hc = 0) is the relaxation time at the critical 
value of the field where the configurations a = +1 and a = —1 are equiprobable (i.e. He = 0). The 
linear response regime breaks down for large ramping speeds when <C Wrev An interesting quantity 
quantifying deviations from the linear-response regime is the fiuctuation-dissipation ratio R defined by, 

R= = , . (52) 

In the limit of small r — > 0, when oc r, then R converges to 1 (in agreement with the fiuctuation- 
dissipation theorem, a result that in the context of steady-state systems has been proven in but 
deviates from 1 as r increases. In the right panel of Figure El we show R[r) when the field is ramped 
from Hi = Q to Hj at different values of Hf. In this case the behavior of both uJ^i^ and R is monotonic 
with r. In Figure El we show the same ramping experiments but comparing, for a given ramping speed, 
the results for t^disj — w^, R as a function of Hf . For values of Hf small enough the ramping process is 
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well described by the linear response-approximation discussed in Sec.jSland 

Let us finish this section emphasizing that the dependence R{r) can be quite complicated and even 
non-monotonic in some cases. Such behavior is observed in the case where the ramping protocol is given 
by H{s) = Ha{^ — '^s/t), i.e. the field starts at a given value Hi = —Ha {Ha denotes the field amplitude) 
and increases until its reversed value Hf = Ha is reached. This case is of much interest regarding heat 
fluctuations and is discussed in detail in Sec. El In Figure [3 we show the behavior of the average work 
W (equal to the average dissipated work as w^ev = due to the independence of the free energy on the 
sign of Ha) and ii as a function of the ramping speed for different values of Ha- 



5 The fluctuation theorem 

Saddle-point equations ()13I14I15|) were derived in the large- limit. Indeed ()2U() is not exact for finite 
but has corrections. However, the results obtained for s{w) exactly satisfy the fiuctuation theorem of 
Crooks 5 . This theorem states the following: let us consider a process where the system is perturbed 
according to a protocol Hp(s) during the time interval [0, i], initially the system being in equilibrium 
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Figure 6: Ramping experiment with Hi = and r = 100 as a function of Hj. The average dissipated work 
(continuous line) and R (dashed line) increase with the field but wjjjg (dotted line) saturates to a finite value equal 
to — log(2) (we represent — u^jjig in order to compare with tf^g). 




Fi gure 7'. Average work w (left) and fluctuation-dissipation ratio R (right) for the case Hj^ — — Hi — Hj^ as a 
function of the ramping speed. The different curves correspond to different values of the amplitude field Ha- These 
are indicated in the plot besides each curve. The straight dashed line in the left panel corresponds to the linear- 
response relation W = 2r in H51|). As explained in last paragraph of Sec. 13 and for this particular protocol, the 
relation = —w'^^ is exact for all values of r and Ha- Moreover, for large values of Ha the work coincides with 
the heat exchanged, see Sec. IHl 
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at the value of the field Hi = Hf{0). We will call this the forward (F) process. Let us now consider 
the reverse process defined as the time-reversed of the forward process: in this process the system starts 
in equilibrium at time s = at the value of the field Hi = Hpit) and the field is changed according to 
the protocol Hr^s) = Hp{t — s). Let the distribution of works generated in this way be Pp(VF), Pr{W) 
for the forward (F) and reverse (R) processes respectively. The two distributions satisfy the following 
relation |5j|, 

where AF = F{Hf) — F{Hi) is the change in the equilibrium free-energy. By rewriting this identity 
as Pr{—W) = Pf{W) exp^ "^^^'^ ^ and integrating it between W = — oo and W = oo leads to the 
Jarzynski equality < exp(— Wdis/ZcsT) >_f= 1 where < ... >f stands for a dynamical average over work 
values obtained along the forward process. 

If we now substitute (|^n|) into the relation we obtain, 

SF{Wd\s) - SR{-Wdis) = (54) 

where we have taken PF,RiW) oc exp(NsF,R{'w)) and we have disregarded the normalization constant in 
the distribution (|^n|) as unimportant. Because the quantity s{w) used in (|^n|) is only exact in the large- 
limit one might be tempted to think that H54|) does not hold. To prove the validity of H54() we rewrite 
(jSSJ in the following way, 

llog{PF{W))-^log{Pn{-W)) = ^^ . (55) 
In the large- limit the distributions ()2U() satisfy, 

lim 1 \og{PF,R{W)) = sf,r{W) (56) 

Af— ►oo iV 

and therefore (|5H) is exact with Wrev = hniAf^oo AF/N. The present approach seems quite general so the 
trajectory entropy derived in a large-A^ theory in any statistical model (interacting or non-interacting) 
must verify the relation 1)54(1 . Another interesting relation that can be obtained from (|54() relates the 
values of for the forward and reverse processes. Differentiating ()54j) respect to w we obtain, 

Mw) + s'r{-w) = \r{-w) - Xf{w) = ^ (57) 

where we used (22). Therefore, the identity s'p{w^^) = Xf{w^^) = (^5]) implies s'^(— it;™P) = 
Xr{—w"^^) = l/ksT. From ()31|) we then infer that for the reverse process is identical to — 
for the forward process and vice versa. This relation is quite interesting because it suggests that in order 
to estimate (e.g. in experiments) the value of w'^ for a given non-equilibrium process it is enough to 
determine tt;™P for the reversed process, a quantity that is experimentally accessible. 

An interesting case of (j5l-{|) occurs whenever the forward and reverse processes are symmetrical mirror 
images, Hf{s) = —Hr{s) = —HF{t — s). This can be accomplished when Ha = Hf = —Hi along the 
forward process and the protocol satisfies H{s) = —H{t — s). In this case the forward and the reverse 
work distributions are identical, W^ev = AF = (or w = Wdis) and (|5l|) reads, 

sH - s{-w) = ^ (58) 

where we have used s{w) = sf{w) = sr{w). The validity of l(55|) can be further demonstrated by close 
inspection of equations (|13|14|15|) . Let s{w) be the value of the dynamical entropy for a given value of the 
work w associated to the value of the Lagrange multiplier A and the magnetization m(s). Then, for the 
reversed value of the work —w, it is possible to show that the corresponding solutions are: —A — I/ZcbT 
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for the Lagrange multiplier and —m{t — s) for the magnetization solution. The resulting entropy is 
then s{—w) = s{w) — w/UbT as given in ((^ . A remarkable consequence of this special case is the 
aforementioned fact that w'^ = — u;™p at any ramping speed and for any value of ther amplitude of the 
field Ha- This case was already shown in Figure [3 The present symmetric case is specially interesting 
because the work done upon the system can be identified with the heat exchanged between the system 
and the bath. The conditions required for such identification are discussed next. 



6 Heat fluctuations and tails 

Until now we focused our efforts to investigate work fluctuations. However, in the same way as the work 
fluctuates, the heat exchanged between the system and the bath also does. The validity of the mechanical 
equivalence of heat (the content of the first law of thermodynamics) suggests that there should not be 
an important difference between heat and work. Heat is more difficult to experimentally measure than 
work and this is the reason why we tend to be more interested in the former. 

A motivation to investigate heat fluctuations has recently arisen in the context of steady state and 
aging systems. In the first case, heat fluctuations were investigated for the simple model of a bead 
dragged through a viscous fluid ^Tj- the second case these were studied for the case of a spin-glass 
model characterized by slow dynamics and aging |18 | I19 | [2(1]. In both cases, a Gaussian component was 
identified in the heat distributions together with some exponential tails. For the steady state system these 
exponential tails were consequence of the validity of an asymptotic fluctuation-theorem for the heat while 
in the aging system the tails were the signature of intermittency effects that have been experimentally 
observed in glasses and colloids |22 [22| ■ 

The heat along a given trajectory can be inferred using energy conservation, —Q + W = AE. To 
extract the heat we just need to know the change in energy AE between the final and initial configurations 
as well as the work W. Here we adopt the sign convention (contrary to that adopted in the introduction 
Sec. that positive heat corresponds to net heat delivered by the system on the surroundings. A 
particular case where work and heat fluctuations are identical is the case described in the preceding 
section where Hj = —Hi = Ha- Due to time reversal symmetry W^cv = = 0. Now, if the field 
amplitude Ha is large enough, the difference in energy AE = —iJ,A{MH) is practically zero so Q = W. 
For example, if Ha = 5, then tanh(5) = 0.99991 (as always we take P = n = 1) so the initial equilibrium 
magnetization is Meq{HA) = —N. The final magnetization after ramping the field to Ha is again of 
order and therefore the fluctuations from trajectory to trajectory in AE are negligible as compared 
to the total work. In Figure |H1 we show the trajectory entropy and free-energy for the case Ha = 10. 
We have chosen to represent variables in terms of heat per particle q = Q/N rather than work to give a 
view of what general shape we can expect from heat distributions. In terms of the heat we expect that 
the same mathematical definitions and relations that we defined in the case of work are also valid. For 
instance, the heat entropy and the heat free-energy are defined in the same way as we did for their work 
counterparts just replacing w by q, in particular J^{q) = q — kBTs{q). Also the equivalent of holds. 

The most probable heat g™P {X{q™^) = 0) and the quantity q"^ (A(g™P) = —1) can also be defined. 
Moreover, a relation equivalent to (|58j) is also expected to hold for large enough values of Ha, 

An interesting aspect of the heat entropy s{q) shown in the left panel of Figure |H1 is the presence of 
quadratic behavior for small values of (; (g ~ 0) together with a linear behavior in the tails {\q\ » 1). 
These characteristic features of the heat entropy s{q) can be inferred by looking at A(g), shown in 
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Figure 8: Trajectory entropy s{q) and trajectory free-energy J^{q) for the case Hf = —Hi = 10 at ramping speeds 
r = 0.1,0.5, 1, 10, 100 (from the most narrower to the most wider distributions). The dashed line in the left panel is 
2/(9) = q/^bT (we take ksT = 1) and is tangent to s{q) at a value . The dashed line in the right panel corresponds 
to y{q) = q and is tangent to T{q) at the value (7™^. Both (7™^, q^ depend on the ramping speed. 

Figure El That figure shows that X{q) is linear with q for q — 0, giving a quadratic behavior for s{q) at 
small values of q. This linear shape in X{q) corresponds to a Gaussian distribution for P{q) = exp(s(g)). 
It also shows that for a wide range of \q\ values there are two plateaus at A(g) ~ A+, — A_ for positive and 
negative values of q respectively. These plateaus correspond to the exponential tails in the distribution. 
This behavior is quite generic at all ramping speeds, the distinction in X{q) between both plateaus and 
the linear behavior at small q becomes more clear as the speed decreases. In such conditions, is not 
very large and the linear response approximation holds. The Gaussian sector describes the statistics of 
small and most probable fluctuations, the exponential tails describe rare events and large deviations. In 
what follows we analyze the Gaussian and exponential tails in more detail. 
In the region where both q, g™P are not too large we have, 

s{q) = -J^^i^ - 9"^)' 9' « 1 (61) 
Substituting this relation in (jSU)) we get, 

implying 

Iq^'cksT ^ ' 

This result shows that the fluctuation-dissipation ratio (|52|) is equal to 1 if heat fluctuations are restricted 
to the sector of q small. Small fluctuations are a key assumption of linear-response theory which also 
leads to (O) . 

This quadratic behavior then goes over straight lines in the most negative and positive sectors of g, 

s{q) = C-\^q q»\ (64) 
s{q)=C^\-q q«-\ (65) 
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Figure 9: X{q) for the same parameters as in Figure|Slin the case r = 0.1. We note the presence of a linear behavior 
for X{q) for small values of q, X{q) = {2/Nal){q - q'^P) {q'^P = 0.207, Na^ = 0.83) and two plateaus for q » 1 
and q « —1 (A+ = 4.95, A_ = 5.95). The former gives rise to the Gaussian component in the heat distribution 
describing the statistics of most probable values. The latter gives rise to two exponential tails for the distribution 
describing the statistics of rare events. 



where C is a constant and A+,A_ correspond to the values of X{q) in the plateaus shown in Figure IHI 
Note that the constant C in ()64I65|) is the same in both sectors. In fact, the relation H60|) imposes such 
constraint between the positive and negative tails in the probability distributions. Substituting 1)64165^ 
into (jEOl) we obtain 

A--A, = J^ (66) 

meaning that the width of the tails is larger for the negative values of q than for positive values. This 
can be interpreted by saying that, despite of the fact that the average heat q is positive, rare fluctuation 
events occur as often for g < (i.e. when the system adsorbs heat from the surroundings) as they 
do for q > (when the system delivers heat to the surroundings). The shape of the heat distribution 
P{q) = exp{s{q)) is then dominated by a central Gaussian distribution with exponential tails at its 
extremes. These features are illustrated in Figure IIUI If the amplitude field is not large enough 
there may be contributions to the heat distribution coming out from the fluctuations in the difference 
in energy between the initial and final configurations. The effect of the value of on the value of the 
average work and the fluctuation-dissipation ratio have been already shown in Figure [71 in particular 
non-monotonic behavior is observed for R. 



7 Finite-size effects 

The method we developed in this paper allowed us to calculate Pn{w) in the large- limit. However, due 
to the non-interacting character of the model, all cumulants of the distribution obtained in the large- 
limit are also exact for finite A^. The proof is quite straightforward. Let us define the generation function 
of all cumulants of the distribution P/v(lF) in (jlj, 

g^(x) = log(exp(xlF)) = log( f dW exp{xW)PN{W)^ . (67) 
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Figure 10: Heat entropy s{q) for the case Hj = —Hi = 10 at ramping speed r = 0.1. Main figure: The sector 
of small or most probable fluctuations g ~ can be well fitted to the Gaussian (|HT|l (dashed line) with parameters 
^mp _ 0.207, A^o"^ = 0.83 satisfying (|63() . The tails extend beyond the Gaussian central part and are of exponential 
type as described in (|64I65|) with A_ = 5.95, A+ = 4.95, C = 9.16. These exponential tails describe the statistics of 
large deviations or rare events. Inset: Heat-distribution P{q) = exp(s(g)) (dots) and the Gaussian fit (the dashed 
line of the main figure) showing that the small q sector of fluctuations (those that are frequently observable) is very 
well fitted by a Gaussian despite of the fact that rare-event tails are big and observable only when plotting s(q) or 
the distribution P{q) in logarithmic scale. 

Cumulants of Pn{W) are obtained using the following formula, 

k being a positive integer. Using the non-interacting character of the model then we can write, 

N 

W = Y,Wi ^ QNix) = Ngi{x) ^ CN{k) = Nci{k) (69) 

i=l 

and therefore all cumulants of the distribution are independent of the size of the system (except by a 
multiplicative constant equal to A'^) . This implies that the expression given for in H48|) and R in 1)52(1 
are independent of N . Therefore, the results we obtained in the large-A^ limit are exact for any finite 
value of N. 

However, albeit cumulants do not depend on N, the shape of the distribution Pn{w) in depends 
on the size N and only in the large- limit the approximate distribution H20|) becomes exact. For 
instance, the value of w"^^ depends on and converges to vJ^ for large enough values of A^. In practice, 
already for N = 5 — 10 convergence of the approximate distribution (j2()jl to the exact result is excellent. 
In order to compare the approximate distributions we obtain from our theory with the exact ones at 
finite A^ we have done numerical simulations of the model. The simulation procedure is described below 
in Sec. |H1 In Figure ^2 we show the distributions we have obtained for A^ = 10 compared to the 
numerical simulations at different ramping speeds. The agreement between theory (continuous lines) 
and simulations (symbols) is good although it worsens progressively as the ramping speed increases and 
the system is strongly driven out of equilibrium. An important feature of the distributions is observed for 
large r and small A^: the presence of a finite fraction of trajectories that dissipate a maximum amount 
of work equal to ffmax = /^(-^/ ~ Hi). For these trajectories the ramping speed is so high and the 
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Figure 11: Distributions Pn{w) for case Hf = —Hi = 10 and = 10 at different ramping speeds (indicated along 
each curve). The continuous hues are the results obtained from the present theory using (|2()|) . The symbols are 
results obtained from numerical simulations of the model for 10^ trajectories. The right panel is the same figure but 
in logarithmic scale. For the largest ramping speed r = 100 there is a finite fraction of trajectories (about 37% of 
the total number of trajectories) where the spins have no time to relax. These trajectories contribute with a singular 
term at it; = li^max = 20 to the distribution PAr=io(u') as described in (f70|) . It cannot be captured by the present 
large- theory so we did not include it in the histogram obtained from the numerical simulations. 

size so small that no change in the initial configuration occurs along the trajectory. We will call these 
trapped trajectories. The fraction of trapped trajectories contributes with a term 6{w — Wmax) to the 
work distribution, 



where Pi\f{w) is a continuous function and a{N) is a size-dependent constant that decreases with N and 
asymptotically vanishes in the large-A'^ limit. The delta function in ()70() is a small N contribution that 
is not captured by the present large- theory. Nevertheless, it might be analytically derived using the 
approach described below in Sec. 17.11 

In Figure El we show the effect of the size on the distributions at a moderate ramping speed. For 
N = 1 the agreement is not good although the behavior of the left tails is reasonably well reproduced. 
However, already for N = 5 the agreement has improved considerably. We conclude that it is between 

= 1 and N = 5 that finite-size effects are important. In Figure ITSl we confirm this strong small-N 
dependence by plotting the most probable work as obtained from the numerical simulations as a function 
of r for different sizes A^ = 1, 2, 5, 20. 

7.1 Reconstructing Pi{w) from the large- theory 

A crucial aspect of the present model is that it is non-interacting. Therefore, if we were to know the 
work probability distribution for A^ = 1 (a single spin) then we could reconstruct the general distribution 
Pn{w). In fact, let Pn{s) denote the Laplace transform of P^iw), 



Pj^{w) = Pn{w) + a{N)6{w - w, 



max 



(70) 




(71) 



Using the result W = J2i=i we can write. 



Pn{s) = {Pi{s)) 



N 



(72) 
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Figure 12: The same as in Figure ITT] but showing the dependence of Pi\f{w) with A'^ at r = 1. The agreement is not 
good for iV = 1 in the central region of the distribution but is reasonably good for the left tail of the distribution 
(see the right plot in logarithmic scale). The agreement improves noticeably beyond A'^ ~ 5. 
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Figure 13: The same parameters as in Figure [TTl but now showing the dependence of the most probable work value 
tD^P with A^. The different symbols in the points correspond to different sizes as indicated in the legend of the figure. 
The continuous dashed line is w^"^ as derived from the theory in the large- limit (where t<;™P = I^dis) the latter 
being independent of A^). In the linear-response regime, r << 1, data have converged to the theory for all sizes. 
Although finite-size effects are large for A^ = 1,2 and r >> 1, already for N = 5 the simulations have converged 
to the theory at all ramping speeds. Data for N = 2 have been connected by a dotted line to emphasize the sharp 
increase of w"^^ around r ~ 5. This sharp increase originates from the presence of two separated peaks in the work 
distribution which height coincide at a given value of the ramping speed around r ~ 5. 
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allowing us to reconstruct Pi\f{w) from the sole knowledge of Pi{w). Although the analytical computation 
of Piiw) might be possible by using other approaches, throughout this paper we have considered a 
thermodynamic approach where the large-A^ theory has been taken as an approximation to finite N. 
This approach turns out to give exact results for all cumulants of the distribution thereby suggesting 
that the reconstruction of Pi{w) from Pi\[{w) might be possible. One could naively think that this is 
possible just using H72() together with the knowledge of Pn{w). Unfortunately, this is not the case as 
the knowledge of PNiw) is only approximate as we showed in the previous section. There is however a 
possible strategy to reconstruct Piiw) that is based on the fact that cumulants are exactly known. Let 
us define the following function, 

hi.) = hm (73) 

N^oo iV 

where gNi^) was defined in (j67|) . In the large- limit we can solve h{x) by applying the saddle-point 
approximation, 

1 



h{x) = lim — log! exp(xVF) 



lim —log([ dW ex.p{xW) exp{Ns{w))') = xw{x) + s{w{x)) . (74) 



N- 

where w{x) is the solution of the equation, 

ds{w 



dw 



= -x . (75) 

w=w(x) 



For a given value of w, (|22|) shows that x = X{w). For instance, for x = 0, —l/ksT we get w = w™^, w'^ 
respectively. Therefore, we can express h{x) in terms of w rather than x, 

h{w) = wX{w) + s{'w) (76) 

By inserting (jMj) in ((75)) we get gi{w) = h(w) and therefore, 

dw' eyi.]i[\{w){w' — w)^Pi{w') = exp^s(ii;)^ (77) 

Formally, this integral equation is closed and provides an exact solution for Pi{w) in terms of the 
entropy s{w). Unfortunately we have been unable to solve it in full generality (as detailed knowledge of 
the solution in ()75|) is required). Yet, for Pi{w) it still holds that there are exponential tails identical to 
those we already derived for P^{w) in the large- A limit. To show this result we use ()22j) and rewrite 
(|77|) as follows, 

j dw' exp (-s(w;) - {w' -w))Pi{w') = l . (78) 

Let us now suppose now that A(w) is approximately constant (equal to A) showing a plateau over a given 

dw 

Pi{w') oc exp(s(u;')) = Cexp(W) , (79) 



region of work values. From (|^^ then s{w') ~ s{w) + ^^^{w' — w) and. 



where C is a constant. This shows that the width of the exponential tail for Pi{w) (and, by extension, 
for Pn{w) at any value of A) is equal to A. 



8 The case of magnetic nanoparticles 

In this section we discuss a system where the previous theory could be experimentally tested. We 
focus our attention on thermally activated magnetic nanoparticle systems |2,'Sj . These systems have 
the great advantage that dynamics is invariant under time-reversal of the magnetic field H — > —H. 
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Also many magnetic field cycles can be experimentally realized in micro-SQUID machines allowing to 
experimentally extract the work distribution with good precision. The main experimental limitation to 
observe WF though is the quite large value of the magnetic moment of the nanoparticle. Transition rates 
are described by the Brown-Neel formula, 

Trelax [H] = Tq exp {^^) (80) 

where tq is a microscopic time describing relaxation within a state and B{H) is a field dependent 
barrier. We consider two cases: A) paramagnetic molecular clusters where the energy barrier is nearly 
field independent B{H) = Bq (this case could also describe specific ferro and ferrimagnetic nanoparticles 
where the anisotropy contribution to the zero-field barrier is negligible, for a discussion see B) 
ferromagnetic nanoparticles with axial anisotropy where B{H) depends on the intensity of the external 
field projected on the easy magnetization axis as described by the Stoner-Wohlfarth expression B{H) = 
Bq{1 — jj- )" where He is the field required to suppress the barrier and a is an exponent in the range 
1.5 — 2. Recent experiments have demonstrated how the height of the barrier Bq can be considerably 
reduced by applying a transverse field, making possible to observe magnetization reversible transitions 
(also called telegraph noise measurements) in single Co nanoparticles at low temperatures |251 126j . 

As we already discussed in Sec. El in a magnetic system a time-reversal invariant protocol can be 
accomplished by switching the magnetic field H from —Ha to H^ {Ha being the amplitude of the field), 
the free energy and the rates being an even function of H. Under such conditions work and heat are 
equivalent if Ha induces a magnetization close to its saturation value. From the experimental point of 
view, it is relevant to understand under which conditions large deviations from the most probable work 
are observable. By large deviations we understand work (heat) fiuctuations corresponding to work (heat) 
values around {o^)- A useful quantity that tell us how difficult it is to sample that region of work 
values is the ratio Vt describing the fraction of trajectories that transiently violate the second law, w <Q. 
This fraction is given by the integrated fiuctuation theorem [S] . This is obtained by rewriting ()53() , 

P„(-H.-) = P„(H.>xp(-^!i^)=e.p(-g!t) . (81) 

where we have taken Pp^iW) = PpiW) = Pr{W). Integrating this expression from W = Q\xptoW = oo 
we obtain, 

^ = ^ n^ = exp(-— — )^>o (82) 
J\l [w > U) kbT 

where M{w < 0),M{'W > 0) are the fraction of trajectories for which the total work is negative and 
positive respectively, 

/O roo 
dWP{W) ; Af{w > 0) = / dWP{W) (83) 
-oo Jo 

and the average on the r.h.s of (|82)) is restricted to the subset of trajectories for which w > 0. Quite 
generally, we expect that O is a non- universal function dependent on all cumulants of s{w), yet its 
exponential dependence in assures that, in the regime where TV are observable, Q is approximately 
described by the value of the average total work divided by the bath temperature W/ksT which is 
approximately given by Nw^^ /ksT . 

We choose Glauber rates as these have been experimentally demonstrated to describe very well the 
relaxation of single magnetic moments |27| I28j . These are given by (|35j) where rreiax(-f^) is given by 
(|8())) . We consider ramping experiments [20] where N particles are subject to the action of a field that 
is switched from H = —Ha up to -fT = Ha at a constant speed H. We generate individual trajectories 
according to the Glauber rates by starting from initial configurations with M = M(.(^{Ha) and repeating 
the ramping protocol many times, each time the total work ^ is computed, W = — /i M{H)dH. 
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Figure 14: Application of the theory to magnetic nanoparticles for a case where w^cv = and w q. Main panel: 
()82p as a function of W /ksT for Cases A and B. Number of particles range from = 1 up to 20 for both cases. 
Case A) corresponds to nanoparticles (open symbols) with /x = SOO/is, T = 200K, Ha = 2T, tq = 10~''s, Bq = 2300K 
at ramping speeds r = lmT/s(circles), lOmT/s (triangles up). Case B) corresponds to ferromagnetic particles (filled 
symbols) with /i = T = 40K, i^^ = 238mT, He = 119mT, Bq = 500K, tq = 10"^, a = 1.5 and ramping speeds 

of O.Ol(circles), 0.1(triangles up), l(diamonds)T/s (continuous, dashed and dotted lines respectively). The continuous 
line is the prediction for a Gaussian work distribution (see discussion at the end of Sec. 0)) in the linear-response 
regime cr^ = 2kBTw^ (it! = 1 in H52|) ). Insets: Both are for ferromagnetic nanoparticles (Case B) with the same 
parameters as in the main panel. Inset a) shows hysteresis cycles for N = 100 ferromagnetic nanoparticles at the 
three ramping speeds (larger hysteresis for higher speeds). Inset b) shows dissipated work distributions at ramping 
speeds O.lT/s for N = 1,5,10,20 particles (larger sizes correspond to narrower distributions). 

If i^sw is the field at which the magnetization of a given particle switches for the first time then, for 
a given trajectory, some of the particles will switch state at a value of the field i/gw < 0, while others 
will switch at H^vi > 0. For fast ramping speeds the dynamics is well described by a first-order Markov 
process |30| and the dissipated work for that trajectory will be identical to the value 2Hs^ averaged over 
all particles. In general, for lower ramping speeds, the relation between the dissipated work and the 
value of HsTK is more complicated. To estimate we generate trajectories and evaluate the fraction of 
them with W > and VF < 0. We chose to do numerical simulations rather than applying the large- A 
theory to give a more clear picture about which results can we expect from a finite number of ramping 
experiments (around 10000). In the main panel of Figure El we plot the value of O (|82|1 as obtained for 
different ramping protocols in cases A and B. All points scatter around a generic (but non- universal) 
curve useful to predict in which regime TV are expected to be observable. An important advantage of 
the time-reversal symmetry property H — s- —H of magnetic nanoparticle systems, as compared to other 
systems jZllHl) is the feasibility of performing many ramping cycles in a single experiment making TV 
observable for Q values as low as 10~^. According to Figure ITU TV should be observable for work values 
as large as 20kBT. 
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9 Conclusions. 



Two-state systems provide a simple conceptual framework to analyze work fluctuations (WF) and tran- 
sient violations (TV) of the second law. These non-equilibrium effects are expected to be relevant and 
observable for nanosized objects when the energies involved are several times ksT, ks being the Boltz- 
mann constant and T the temperature of the bath. These have been already observed in the unfolding of 
small RNA hairpins [7] as well as in polysterene beads dragged through a solvent |S]. Related measure- 
ments include the experimental test of the Gallavotti-Cohen fluctuation theorem in Rayleigh-Bernard 
convection ^ and turbulent flows [^. Other experiments include the observation of gravitatory po- 
tential energy fluctuations in driven granular media [2^. The scientific discipline behind all such rich 
phenomenology deserves to be called thermodynamics of small systems. It deals with the thermal be- 
havior of non-equilibrium small systems where the typical energies are few times ksT. The statistics 
of energy exchange processes between the system and the thermal environment is described by frequent 
Gaussian distributed events plus rare events corresponding to large statistical deviations from the average 
value. The theoretical and experimental study of these fluctuations could be of relevance to understand 
issues related to the organization and function of biological matter in the nanoscale [221 • 

In this paper we studied WF in two-state systems. We have introduced a trajectory thermodynamics 
formalism with the specific aim to quantify WF in such model. We have shown how to define a trajectory 
entropy s{w) that characterizes WF around the most probable value w"^^, and a trajectory free-energy 
J-{w) whose minimum value sA w = w'^ specifies the value of the work that needs to be efficiently 
sampled to quantitatively test the Jarzynski equality. The theory requires the introduction of a Lagrange 
multiplier \{w), its inverse playing the role of a temperature in the trajectory thermodynamics formalism. 
Analytical expressions for the trajectory potentials s{w), J-{w) have been also derived. In general, both 
values tD^P and are of the same magnitude but opposite sign, meaning that large deviations of WF 
need to be sampled to recover equilibrium free-energies from non-equilibrium measurements, e.g. by 
using the Jarzynski equality. 

We have then carried out a systematic study of WF in the framework of the large- theory. Several 
results are worth mentioning. First of all, we have found an analytical expression for the trajectory 
entropy that satisfies the fluctuation theorem by Crooks [S] that relates forward and reverse processes. 
An important result is that the value of the work that has to be sampled in order to test the 
Jarzynski equality is equal to the most probable value of the work (with a minus sign) for the reverse 
process. Intuitively this means that the forward and reverse distributions must overlap each other in 
order to get good estimates of the work using the Jarzynski equality, a result that was emphasized 
long-time ago by Bennett j33j . Furthermore, if both forward and reverse processes are symmetric mirror 
images then t/;''^^ = and = —w'^^ independently of how far the system is driven out of equilibrium 
. This last case is particularly interesting as the total work practically coincides with the heat. The 
fluctuation theorem by Crooks is then also applicable to the heat in that limit, a result that is quite 
reminiscent of a heat fluctuation theorem recently derived |17| I34 j . For the heat distribution, we find 
that it is described by a central Gaussian distribution describing local equilibrium, i.e. with R = 1, and 
long exponential tails with widths described by the Lagrange multiplier \{w), which plays the role of 
the inverse of a temperature. Strictly speaking, because the temperature must be a positive quantity, 
only the tails in the negative sector g << — 1 where A is negative admit such an interpretation (i.e. in 
the sector of WF dominated by TV). It has not escaped our attention that this temperature could be 
related to other non-equilibrium temperatures that have been defined in other contexts [351 136j . such as 
steady-state j^Zj or aging systems jTS]. 

Our study raises the following question: to what extent are work and heat fluctuations equivalent? 
We already emphasized in Sec. El that work and heat should be equivalent, at least this is the underlying 
content of the flrst law of thermodynamics. However, from the perspective provided by the present 
analysis, some important differences can be underlined. Exponential tails are more often observed in 



26 



the heat rather than in the work. Such result has been exphcitly shown in the case of a bead dragged 
through a fluid ^7] where the work is clearly Gaussian distributed while the heat displays exponential 
tails. However, in that case the origin of this difference lies on the fact that the motion for the bead is 
described by a stochastic linear equation which in general might not be the case. The difference between 
heat and work has its root at the true microscopic definition of these quantities. Heat is identical to work 
when the final energy of the system is constrained be identical to the initial value (i.e. Q = W ii AE = 0, 
for the heat we adopt the sign convention of Sec. E]). The simplest interpretation is that exponential 
tails in the work distribution are always present if the model is non-linear by definition (which is not the 
case for the aforementioned case of the bead dragged through the fluid). However, work distributions 
always tend to be masked by a Gaussian contribution coming out from the Gaussian fluctuations that 
characterize the initial equilibrium state. Therefore, only when thermal fluctuations in the initial and 
final states are negligible as compared to the total amount of work along the trajectory, the measured 
work distributions are paralleled by the heat distributions and tails can be observed. This explains the 
qualitative difference observed between the functions X{w) in Figs. 1^1 and the right panel in Fig. |3J In 
the latter, Gaussian fluctuations in the energy of the initial and final configurations tend to mask the 
presence of the exponential tails. 

We also studied finite-size effects to test how good the large- theory is and provided a strategy to 
re-derive the finite-A^ work distribution from the large-A^ result. An important conclusion is that the 
large-A^ theory accounts for the existence of exponential tails also at finite A^, the value of the widths 
A+,A_ (corresponding to the plateaus in X{w)) being independent of A^. In addition, we applied the 
theory to magnetic nanoparticle systems which provide an experimental realization of two-state systems. 
We studied under which conditions the theory can be experimentally tested. Our results suggest that 
WF and TV should be observable whenever average work values are not much larger than 20/cbT. It 
is realistic to say that we are currently at the limit of the resolution of current micro-SQUID devices 
for the detection of single small magnetic moments (around few hundreds of ^b)- Surely, we will see 
developments in the near future and experimental measurements of WF in magnetic systems, as well as 
the test of the present theory, will become possible. 
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